94 цифры числа e

Реализация алгоритма нахождения 94 цифр числа e.

Предыстория

Первый вариант реализации вычисления цифр числа e указан в журнале Наука и Жизнь №4 1990 г.: вырезка в формате djvu. Там получилось только 27 цифр. Через год с небольшим опубликовали вариант, где цифр уже стало 92: Наука и жизнь №6 1991.

По сути, это уже максимум. При вычислениях (нахождение остатка от деления), чтобы не потерять точность при переполнении, приходится хранить только 6 цифр, а не 7. В результате получаем 15 регистров × 6 знаков = 90. Начало используется с переполнением, поэтому 91 цифра после запятой или всего 92.

Реализация

Текущая реализация – это немного улучшенный вариант того, что было опубликовано в 1991:

• Программа самостоятельно обнуляет все регистры перед вычислением, что позволяет запустить её многократно, например, для определения производительности ПМК.
• Она ещё и короче на 31 команду.

С учетом того, что в конце вычислений уже можно не боятся переполнений, в последний регистр дописаны ещё две цифры числа e. Так и получается 94 цифры.

Глубина расчёта

Сначала поясню загадочное число 65, используемое в программе, но мало освещённое в исходной статье. Дело в том, что точность суммы сходящегося ряда (суммы), определяется последним членом. При вычислении e по формуле ряда ∑ 1⁄N!, точность определяется последним 1⁄N! Для определения разрядности (десятичной) числа используется десятичный логарифм. Количество разрядов у числа (точность) определяется значением логарифма. Если учесть, что N! = 1×2×3×…×N, а логарифм произведения равен сумме логарифмов отдельных членов, то достаточно вычислить lg(1)+lg(2)+…+lg(N) для получения количества знаков у числа N! Например, если результат больше трех, то число содержит больше трех знаков (четыре).

Для вычислений суммы логарифмов составим простую программу:

На входе число N, на выходе сумма логарифмов. Зададим 4 и получим 1.3802112, т.е. 2 знака, как и ожидалось: 1×2×3×4 = 24.

А если задать 65, то получим 90.916331, те самые 91 знак после запятой. Вот откуда это число 65 – это предел суммы ряда, чтобы получилась 91 точная цифра.

Алгоритм расчёта

Используя разложение в цепную дробь, на каждом шаге вычисляем 1/N. Для дальнейшего понимания, как делим на N частями, используем простой пример. Путь хотим вычислить 1/73. Возьмём вместо единицы 10⁸ и поделим группами по 2 знака (10²=100).

Таким образом, получим цифры результата (из 3-й колонки) 01369863 или с учётом 10⁸ это 1.3698630-02. Если вычислить 1/73 напрямую, то получим то же.

В нашем случае вместо единицы будем использовать 10⁹¹. Делать будем целочисленно группами по 6 цифр, первая группа может содержать 7. Получится 10⁷ и 14 раз по 10⁶. В каждой группе результат (целый, 6 знаков) сохраняется в регистре, а остаток переносится. Таким образом пробегаем по всем группам (регистрам). Итог деления - в регистрах, и остаток в стеке. Повторяем N (65) раз с уменьшением N.

Алгоритм целочисленного деления и нахождения остатка выполняется в подпрограмме, а накопленный остаток хранится в стеке. В исходной программе делается многократный повтор: извлечь из очередного регистра текущие разряды, вызвать подпрограмму, сохранить. После пробега по всем разрядам уменьшить число N. При обнулении – закончить.

Оптимизация

Я немного поменял роли. Остаток по-прежнему в стеке, но регистр R0 стал использовать для хранения числа N. Поэтому вычисляемые знаки, хранимые ранее в R0, тоже сохраняются в стеке в регистре RT. Какой вместительный стек :). Это позволило даже чуть-чуть сократить манипуляции в стеке, т.к. сохранять N уже не нужно, и конец цикла сделался проще.

Кроме того, чтобы сократить длину программы, многократный вызов подпрограммы с разными регистрами сделал во внутреннем цикле. А для счётчика этого цикла использовал два последних разряда регистра Re, потому что для вычислений используется только 6 разрядов.

Микрооптимизация: для зануления 15 регистров вместо явного числа 15 (две команды), используется дробное число 1/65 (функция обратной величины - одна команда), у которого в конце как раз 15. А для косвенной адресации без разницы.

По окончании хранимое в стеке начало числа e (в оригинале R0), нужно сбросить в R0. Да и в регистре Re нужно почистить хвост. Тут я немного сжульничал. Вместо прокрутки стека (из RZ в R0) и прибавления единицы, я вычислил средствами ПМК число e (как раз те первые 7 цифр, что в RT, уже с +1). Желающие могут проверить, что в стеке в RZ, RT то же самое, что в R0 (без +1). И вместо обнуления двух последних разрядов Re дописал два настоящих разряда числа e.

Результат

Как и в исходной программе результат хранится в регистрах R0…Re. Начало (R0) продублировано на экране (RX). Всего 94 цифры.

Для повтора вычислений нажмите В/0 С/П.

Текст программы

Детали программы

Подпрограмма с адреса 52

На входе в стеке предполагается:

На выходе должны получить:

Кстати, небольшие уточнения. Первое, на входе RX и RZ можно поменять местами, что и делается в самом начале (первый вызов подпрограммы). Второе, делитель состоит из 2 знаков, остаток может тоже. При умножении на 10⁶ уже 8 знаков. Это предел для ПМК. Поэтому группы по 6 знаков, а не по 7. В первой группе явно начинаем с 10⁷, поэтому её это ограничение не касается. В первой группе может быть, а в конце и будет, 7 знаков.

Теперь по порядку, что выполняем в подпрограмме:

1. Удаляем M, путём откидывания в хвост стека. Адреса 52..53.
2. K+S продвигаем в стек и вычисляем [(K+S)/N]. Адреса 54..55.
3. Из продвинутого выше K+S вычитаем N×[(K+S)/N], получаем остаток. Адреса 56..61.
4. Результат проверяем на ноль. Нулём будет, как минимум, в самом конце, потому что при N=1 делится без остатка. Заодно проверка делает результат устойчивым к операции ВП6 (см. X2-влияющая команда).
5. Благодаря тому, что команда ВП регистр X1 не трогает, где осталось N×[(K+S)/N], мы это вытаскиваем (адрес 66), делим на N (точно поделится) и получаем то, что нужно.

Начало программы

N (=65) загоняем в R0. Теперь нужно очистить 15 регистров. Это делают команды по адресам 06..10, но в R1 будет не 15, а 1/65 = 1.5384615^-02, которое и содержит 15 (в конце). Для косвенной адресации это не важно (см. Число положительное, но меньше единицы). Цикл очистки прервется, когда занулится сам R1.

Общий (внешний) цикл

С адреса 11. Как упоминалось, Re используется как счетчик внутреннего цикла по всем регистрам. Тут предполагается, что RX = Re×100. В самом начале это ноль, что верно.

Начальный индекс 1 сохраняется в Re. Чтобы не потерять из RT число H, 10⁷ получаем через ВП без сдвига стека. После Cx это безопасно. Напомню, в RZ хранится первая группа H (в самом начале ноль), RY забили мусором (новое Re), и начинаем с 10⁷. Это как раз случай, где входные RX, RZ для подпрограммы идут в другом порядке.

После вычисления нового H командами по адресам 19..20 загоняем его в RZ. Число S остается в RX и RY. Дубль S потом потеряется.

Теперь внутренний цикл по всем регистрам начиная с адреса 21. Извлекаем очередные 6 знаков, вызываем подпрограмму, сохраняем (адреса 21..24). Теперь нужно увеличить Re на 1. Чтобы не портить стек, превращаем число (всегда положительное) в 1 через КЗН. Новое Re сохраняем (команды 25..28).

Чтобы определить, что внутренний цикл закончен, косвенно извлекаем по Re. Если это само Re и есть, то дошли до конца. Что и делает проверка по адресу 31.

Для Re нужно отрезать две последние цифры, что делается по адресам 34..37. Здесь и далее сознательно используем опасную в программном режиме ВП, вместо F10^x, чтобы не сдвигать стек и не потерять значение H. После П→xe команда ВП безопасна.

После последнего вызова подпрограммы умножаем на 100 новое Re. Здесь, после возврата из подпрограммы, команда ВП опять безопасна. И возвращаемся на начало внешнего цикла, если не кончился N.

Окончание

При выходе из цикла (адрес 44), с учётом того, что Re уже умножен на 100, последние цифры числа e прибавляем. И делаем жульничество, описанное в разделе по оптимизации.

Всё.